function [] = main_Server(Tn, nobs_in, qn_in)

% Main simulation driver for Monte Carlo experiments.
% Depending on the setting of interest, users can uncomment different parts of the code:
% - Simpler examples with constant risk premia (as in Figure 1 of the main text)
% - No jump risk components (as in Figure 1 of the main text)
% - Full model with jumps and time-varying premia (as in the appendix)

% Suggested configurations:
% Tn = 5, 10, 15, 20, 25, 30, 35, 40, 45, 50;
% Corresponding number of observations per block:
% nobs = [78, 26, 13, 26/3, 13/3, 13/7, 1, 1/3, 1/7, 1/21];
% Corresponding values for qn_in = [21, 42, 63, 126];

% Example usage:
% main_Server(5, 78*21, 21);   



rng(100)

% This code illustrates how the RMSE varies with Tn (total time), nobs (number of observations per day), and qn (window size).
% Note: `nobs_in` is the number of total intraday observations per asset; 
nobs = nobs_in / 21;

MC = 1000;            % Number of Monte Carlo replications
jump = false;         % Whether to include jump detection in estimation. 'true': used for debugging purpose.
K = 1;                % Number of factors
l = 1;                % 

qn = round(nobs * qn_in);     % Window length for estimation
Delta_n = 1 / 252 / nobs;     % Observation interval at current sampling frequency
Delta_n_DGP = 1 / 252 / 78;   % DGP is generated at 5-minute frequency (78 obs/day)

% Constraint: This simulation assumes the DGP is always generated at 5-minute frequency,
% then aggregated to match other sampling frequencies. 
% Hence, the current code only works when nobs ≤ 78.
% To generalize, replace '78' with the desired base frequency.

n_DGP = floor(Tn / Delta_n_DGP);   % Number of DGP steps
n = floor(Tn / Delta_n);           % Number of observed data points (after aggregation)

M = 25;            % Number of assets (cross-sectional dimension)
a = 0.02;          %
b = 0.5;           %

L_k = ones(1, K);  %
H = sum(L_k);
p_g1 = 1.0;%0.7; % L_m = [0,T]
p_g2 = 0;%0.1; % L_m = [\zeta_m, T]
p_g3 = 0;%0.1; % L_m = [0, \theta_m]
p_g4 = 0;%0.1; % L_m = [\zeta_m, \theta_m]
M_g1 = floor(M*p_g1);
M_g2 = floor(M*p_g2);
M_g3 = floor(M*p_g3);
M_g4 = floor(M*p_g4);


% Tuning Parameters
un = 3;
vn = 1000*log(1/(Tn*Delta_n));

% Parameters for the variance process

sigma_I = rand(M,1)*0.1 + 0.05;
kappa_sigma = 3;
theta_sigma = 0.06*ones(K,1);
sigmaF_0 = theta_sigma;
gamma_sigma = 0.3*ones(K,1);
rho = -0.5; % correlation of factor and volatility process (leverage effect)
rho_ls = 0.5; % correlation between risk premia and sigma process

% Parameters for the risk premia process
%lambda_C  =  randn(K,1)*0.01; % 
%lambda_J = randn(H,1)*0.01; 
% Continuous risk premia
kappa_lambda_C = ones(K,1);
theta_lambda_C = 0.1;%0.05*(rand(K,1)-0.5)*2;

lambda_C  =  theta_lambda_C; %initial value of continuous lambda

gamma_lambda_C = 0.1;
% Jump risk premia
kappa_lambda_J = ones(H,1);
theta_lambda_J = zeros(H,1)-0.05*4; 
% for i = 1:K
%     theta_lambda_J(2*i-1,:) = 0.05;%*rand(1,1); % negative jump risk premia
%     theta_lambda_J(2*i,:) = (-0.05);%*rand(1,1);
% end

lambda_J = theta_lambda_J; % initial value of jump lambda

gamma_lambda_J = 0.1;

% Parameter for the beta process
beta1_C    = 0.5*randn(M,K)+[ones(M,1), zeros(M,K-1)];
% continuous beta
kappa_beta_C = 2;
theta_beta_C = [ones(M,1), zeros(M,K-1)];
%theta_beta_C = [ones(M,1)+randn(M,1), randn(M,K-1)];
gamma_beta_C = 0.5;
% jump beta: stack all jump betas together
beta_J = b*(rand(M,H)-0.5)*2+[ones(M,L_k(1)), zeros(M,H-L_k(1))];

% Parameters for the jump process
jp_trunc_F = 0.1; % truncation thresholds for the factors
jp_trunc_P = 0.2; % truncation thresholds for the individual stocks

B_0  = repmat([0, 0], K, 1);
B_kl = [-inf, inf];

nu_F = 63; 
nu_I = randi([floor(nu_F*2/3),floor(nu_F*4/3)],M,1);


p_F = 0.5*ones(K,1);%
p_I = rand(M,1);%positive jump probability

JSize_F = 0.01;

JSize_P = 0.1*sqrt(sigma_I);

Avg_jp = zeros(H,1);
pd_F = makedist('Exponential', JSize_F); % create exponential distribution
t_F = truncate(pd_F, 0, jp_trunc_F); % create truncated exponential distribution
for i = 1:K
    for j = 1:L_k(i)
        m = sum(L_k(1:i-1));
        fun = @(x) pdf(t_F,x).*x;
        if min(B_kl(m+j,:)) >= 0
            Avg_jp(m+j,1) = p_F(i,1)*integral(fun, min(B_kl(m+j,:)), max(B_kl(m+j,:)));
        elseif max(B_kl(m+j,:)) <= 0
            Avg_jp(m+j,1) = -(1-p_F(i,1))*integral(fun, -max(B_kl(m+j,:)), -min(B_kl(m+j,:)));
        else
            Avg_jp(m+j,1) = p_F(i,1)*integral(fun, 0, max(B_kl(m+j,:))) - (1-p_F(i,1))*integral(fun, 0, -min(B_kl(m+j,:)));
        end
    end
end


Avg_B0_jp = zeros(K,1);

for i=1:K
   
    Avg_B0_jp(i, 1) =  p_F(i,1)*integral(fun, 0, max(B_0(i,:))) - (1-p_F(i,1))*integral(fun, 0, -min(B_0(i,:)));
    
end


% Estimator Initialization
LambdaT   = zeros(K+H,MC); % average risk premium
LambdaT_C = zeros(K,MC);
LambdaT_J = zeros(H,MC);
Lambda_hat_AJX = zeros(K+H,MC); % AJX estimate of Lambda_T
Vn_AJX = zeros(K+H, K+H, MC);
Lambda_hat = zeros(K,MC); % AJX estimate of Lambda_T
Vn_hat = zeros(K, K, MC);
Lambda_hat_Avg  = zeros(K,MC);
beta_C_hat_AJX = nan(M,K,(floor(n/qn)-1),MC);
beta_C_tr = nan(M,K,(floor(n/qn)-1),MC);

beta_J_hat_AJX = nan(M,H,MC);
Lambda_hat_AJX_truebeta = zeros(K+H,MC);

tic;
for iMC = 1:MC
    rng((l-1)*MC+iMC)
    disp(iMC);
    % Simulation for the variance process
    % Generate the Brownian motion

    dW_sigma = randn(K,n_DGP);
    dW_F = sqrt(1-rho^2)*randn(K,n_DGP)+rho*dW_sigma;
    % Initialize the matrix
    sigmaF_t = zeros(K,n_DGP);
    sigmaF_t(:,1) = sigmaF_0;
    dFt_inov = zeros(K,n_DGP);
    % Setting initial value
    dFt_inov(:,1) = sigmaF_0.^0.5.*randn(K,1)*sqrt(Delta_n_DGP);
    % CIR process for variance
    for i = 2:n_DGP
        sigmaF_t(:,i) = abs(sigmaF_t(:, i-1) + kappa_sigma*(theta_sigma-sigmaF_t(:,i-1))*Delta_n_DGP ...
                        + gamma_sigma.*sqrt(sigmaF_t(:,i-1)).*dW_sigma(:,i)*sqrt(Delta_n_DGP));
        dFt_inov(:,i) = sigmaF_t(:,i-1).^0.5.*dW_F(:,i)*sqrt(Delta_n_DGP);
    end
    
    % Simulation for the Jump part
    PoissonF = binornd(1,nu_F*Delta_n_DGP*ones(K,n_DGP)); % random jumping vector
    for i = 1:K
        ix_bad = PoissonF(i, 1:end-1)==1 & PoissonF(i, 2:end)==1;
        PoissonF(i,ix_bad) = 0;
    end
    
    JFsgn = binornd(1,p_F.*ones(K,n_DGP));
    JumpVectorF = (JFsgn.*random(t_F, K, n_DGP) + (1 - JFsgn).*(-random(t_F, K, n_DGP))) .* PoissonF;
    JumpVectorF_compen = JumpVectorF - Delta_n_DGP*nu_F*(p_F*mean(t_F) - (1-p_F)*mean(t_F));
    
    PoissonP = binornd(1,Delta_n_DGP*nu_I.*ones(M,n_DGP));
    for i = 1:M
        ix_bad = PoissonP(i,1:end-1) == 1 & PoissonP(i,2:end) == 1;
        PoissonP(i,ix_bad) = 0;
    end
    JPsgn = binornd(1,p_I.*ones(M,n_DGP));
    JumpVectorX = zeros(M,n_DGP);
    JumpVectorX_compen = zeros(M,n_DGP);
    for i = 1:M
        pd_P = makedist('Exponential', JSize_P(i));
        t_P = truncate(pd_P, 0, jp_trunc_P);
        JumpVectorX(i,:) = (JPsgn(i,:).*random(t_P, 1, n_DGP) + (1-JPsgn(i,:)).*(-random(t_P, 1, n_DGP))).*PoissonP(i,:);
        JumpVectorX_compen(i,:) = JumpVectorX(i,:) - Delta_n_DGP*nu_I(i).*(p_I(i)*mean(t_P) - (1-p_I(i))*mean(t_P));
    end
    
    % Simulation for the continuous risk premia and beta
    beta_C    =  zeros(M,K,n_DGP);
    beta_C(:,:,1) = beta1_C;
    lambdat_C  = zeros(K,n_DGP); 
    lambdat_C(:,1) = lambda_C;
    lambdat_J = zeros(H,n_DGP);
    lambdat_J(:,1) = lambda_J;
    
    
    rns = randn(M*K+K+H,n_DGP);    
    rns((M*K+1):(M*K+K),:) = sqrt(1-rho_ls^2)*rns((M*K+1):(M*K+K),:) + rho_ls*dW_sigma; % correlation between lambda and c^F    
    rns_ = nan(M*K+K+H,n_DGP);
    rns_((M*K+1):(M*K+K+H),:) = rns((M*K+1):(M*K+K+H),:);
        
    rho_bl = rand(M,K); %correlation between beta and lambda
        
    for j = 1:K
       
        rns_((M*(j-1)+1): M*j,:) =  sqrt(1- repmat(rho_bl(:,j),1,n_DGP).^2) .* rns((M*(j-1)+1): M*j,:)+  repmat(rho_bl(:,j),1,n_DGP).* repmat(rns(M*K+j,:),M,1);  % Mx1 M x n
                
    end
    rns = rns_;
    
    for j = 1:K
        for i = 1:M
            u = 1-kappa_beta_C*Delta_n_DGP;
            v = gamma_beta_C*sqrt(Delta_n_DGP);
            w = kappa_beta_C*theta_beta_C(i,j)*Delta_n_DGP;
            beta_C(i,j,:) = filter(1, [1 -u], v*rns((j-1)*M+i,:)+w, beta_C(i,j,1));
            %beta_C(i,j,:)= repmat(beta_C(i,j,1),1,n); %constant beta
        end
        u = 1 - kappa_lambda_C(j)*Delta_n_DGP;
        v = gamma_lambda_C*sqrt(Delta_n_DGP);
        w = kappa_lambda_C(j)*theta_lambda_C(j)*Delta_n_DGP;
        lambdat_C(j,:) = filter(1, [1 -u], v*rns((M*K+j),:)+w, lambda_C(j));
       % lambdat_C(j,:)  =  repmat(lambda_C(j),1,n); % constatn risk premia
    end
    
    for j = 1:H
        u = 1 - kappa_lambda_J(j)*Delta_n_DGP;
        v = gamma_lambda_J*sqrt(Delta_n_DGP);
        w = kappa_lambda_J(j)*theta_lambda_J(j)*Delta_n_DGP;
        lambdat_J(j,:) = filter(1, [1 -u], v*rns((M*K+K+j),:)+w, lambda_J(j));
        %lambdat_J(j,:) = repmat(lambda_J(j),1,n);% constatn risk premia
    end
    
    LambdaT_C(:,iMC) = mean(lambdat_C,2);
    LambdaT_J(:,iMC) = mean(lambdat_J,2);
    LambdaT(:,iMC) = [LambdaT_C(:,iMC); LambdaT_J(:,iMC)];
    
    % Generate factor dynamics
    dFt_tilde = zeros(K,n_DGP);
    dFt_bar = zeros(H,n_DGP);
    for j = 1:K
        %u = JumpVectorF(j,:) > min(B_0(j,:)) &...
        %    JumpVectorF(j,:) < max(B_0(j,:));
        dFt_tilde(j,2:n_DGP) = lambdat_C(j,1:(n_DGP-1))*Delta_n_DGP;% + JumpVectorF(j,:).*u;
        for k = 1:L_k(j)
            m = sum(L_k(1:j-1));
            dFt_tilde(j,2:n_DGP) = dFt_tilde(j,2:n_DGP) + lambdat_J(m+k,1:(n_DGP-1))*Delta_n_DGP;
        end
    end
    dFt_tilde = dFt_tilde + dFt_inov;
    dFt = dFt_tilde;
    for j = 1:K
        for k = 1:L_k(j)
            m = sum(L_k(1:j-1));
            u = JumpVectorF(j,:)> min(B_kl(m+k,:)) &...
                JumpVectorF(j,:)< max(B_kl(m+k,:));
            dFt_bar(m+k,:) = JumpVectorF(j,:).*u;
        end
    end
    dFt = dFt+ JumpVectorF_compen;
    
    % Simulate Stock Prices
    dPI_t    = diag(sigma_I)^0.5*randn(M,n_DGP)*sqrt(Delta_n_DGP) + JumpVectorX_compen;          % residual innovations: corresponds to sigma^I dW
    dPt    = zeros(M,n_DGP); % return on asset
    
    u = zeros(K,n_DGP);
    for j = 1:K
        u(j,:) = JumpVectorF(j,:)>min(B_0(j,:)) & ...
                 JumpVectorF(j,:)<max(B_0(j,:));
    end
    
    for t = 2:n_DGP
        dPt(:,t) = beta_C(:,:,t-1)*(lambdat_C(:,t-1)*Delta_n_DGP + dFt_inov(:,t) + (JumpVectorF(:,t).*u(:,t)- Delta_n_DGP*nu_F*(Avg_B0_jp)));
    end
    dPt = dPt + beta_J*lambdat_J*Delta_n_DGP + beta_J*(dFt_bar  - Delta_n_DGP*nu_F*repmat(Avg_jp,1,n_DGP)) + dPI_t;
    
    % Generate incomplete data
    zeta_mat_g2 = randi([1,n_DGP], M_g2, 1);
    filter_mat_g2 = NaN(M_g2, n_DGP);
    for i = 1:M_g2
        filter_mat_g2(i, zeta_mat_g2(i):end) = 1;
    end
    
    theta_mat_g3 = randi([1,n_DGP], M_g3, 1);
    filter_mat_g3 = NaN(M_g3, n_DGP);
    for i = 1:M_g3
        filter_mat_g3(i, 1:theta_mat_g3(i)) = 1;
    end
    
    rand_mat_g4 = randi([1,n_DGP], M_g4, 2);
    zeta_mat_g4 = min(rand_mat_g4, [], 2);
    theta_mat_g4 = max(rand_mat_g4, [], 2);
    filter_mat_g4 = NaN(M_g4, n_DGP);
    for i = 1:M_g4
        filter_mat_g4(i, zeta_mat_g4(i):theta_mat_g4(i)) = 1;
    end
    
    filter_mat = [ones(M_g1,n_DGP); filter_mat_g2; filter_mat_g3; filter_mat_g4];
    
    dPt = dPt.*filter_mat;
    cdPt = cumsum([zeros(M,1),dPt]')';
    dPt = diff(cdPt(:,1:round(n_DGP/n):end)')';
    cdFt = cumsum([zeros(K,1),dFt]')';
    dFt = diff(cdFt(:,1:round(n_DGP/n):end)')';

    cJumpVectorF = cumsum([zeros(K,1),JumpVectorF]')';
    JumpVectorF =  diff(cJumpVectorF(:,1:round(n_DGP/n):end)')';%JumpVectorF(:,1:round(n_DGP/n):end);
    cJumpVectorX = cumsum([zeros(M,1),JumpVectorX]')';
    JumpVectorX = diff(cJumpVectorX(:,1:round(n_DGP/n):end)')';
    
    beta_C = beta_C(:,:,1:round(n_DGP/n):end);
    % AJX Estimator
    [Lambda_hat(:,iMC), Vn_hat(:,:,iMC)] = FM(dPt, dFt, Delta_n);
    [Lambda_hat_AJX(:,iMC), Vn_AJX(:,:,iMC), beta_J_hat, beta_C_hat,Lambda_hat_AJX_truebeta(:,iMC)]  = AJX_HF(dPt, dFt, L_k, B_kl, Delta_n, un, qn, vn, JumpVectorF, JumpVectorX, jump,beta_C,beta_J);

    beta_C_hat_AJX(:,:,:,iMC) = beta_C_hat;
    beta_J_hat_AJX(:,:,iMC) = beta_J_hat;
    Lambda_hat_Avg(:,iMC) = mean(dFt,2);
    beta_C_tr(:,:,:,iMC) = beta_C(:,:,1:qn:end-qn);
   
end
toc;

%         



save(['Simu_data_Tn_',num2str(Tn),'_nobs_',num2str(nobs_in),'_qn_',num2str(qn_in),'.mat'],'LambdaT_C','LambdaT_J','Lambda_hat','Lambda_hat_AJX', 'Vn_AJX',  'Vn_hat', 'LambdaT', 'MC', 'Tn', 'nobs', 'Delta_n', 'n',...
    'M', 'K', 'L_k', 'B_kl', 'H', 'un', 'qn', 'vn', 'sigmaF_0', 'sigma_I', 'kappa_sigma', 'theta_sigma', 'gamma_sigma',...
    'rho', 'lambda_C', 'lambda_J', 'kappa_lambda_C', 'theta_lambda_C', 'gamma_lambda_C', 'kappa_lambda_J', 'theta_lambda_J',...
    'gamma_lambda_J', 'beta1_C', 'kappa_beta_C', 'theta_beta_C', 'gamma_beta_C', 'beta_J', 'nu_F', 'nu_I', 'jump', 'beta_J',...
    'p_g1', 'p_g2', 'p_g3', 'p_g4', 'a', 'b','beta_J_hat', 'Avg_jp','Avg_B0_jp','Lambda_hat_Avg','beta_C_hat_AJX','beta_J_hat_AJX','beta_C_tr','Lambda_hat_AJX_truebeta');

end

